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Abstract 

Sea-ice  motion  consists  of  complex  non-rigid  motions  in¬ 
volving  continuous,  piece-wise  continuous  and  discrete  par¬ 
ticle  motion.  Techniques  for  estimating  non-rigid  motion 
of  sea  ice  from  pairs  of  satellite  images  (generally  spaced 
three  days  apart)  are  still  in  the  developmental  stages.  For 
interior  Arctic  and  Antarctic  pack  ice,  the  continuum  as¬ 
sumption  begins  to  fail  below  the  5  km  scale  with  evidence 
of  discontinuities  already  revealed  in  models  and  remote 
sensing  products  in  the  form  of  abrupt  changes  in  magni¬ 
tude  and  direction  of  the  differential  velocity.  Using  a  hi¬ 
erarchical  multi-scale  phase-correlation  method  and  prof¬ 
iting  from  known  limitations  of  cross  correlation  methods, 
we  incorporate  the  identification  of  discontinuities  into  our 
motion  estimation  algorithm,  thereby  descending  below  the 
continuum  threshold  to  examine  the  phenomenon  of  discon¬ 
tinuous  non-rigid  sea-ice  motion. 

1.  Introduction 

Sea-ice  consists  of  a  collection  of  floes  1  with  differential 
motion  subject  to  non-local  forcing  at  an  aggregate-scale 
with  ranges  from  10-300  km  spatially  and  1-20  days  tem¬ 
porally  [37],  [34],  This  differential  motion  results  in  fea¬ 
tures  such  as  leads,  families  of  slip  lines,  cracks,  and  ridges. 
Within  the  sea-ice  community  there  is  no  formal  definition 
of  this  scale.  Coon,  et  al.  [14],  Aksenov  and  Hibler  [1], 

1 A  single  piece,  large  or  small,  of  ice  formed  on  the  surface  of  an  ocean 
or  sea. 


and  Wilchinsky  and  Feltham  [42]  have  made  progress  in 
developing  theoretical  “anisotropic”  sea-ice  rheologies,  but 
adequate  data  products  to  rigorously  test  their  theories  do 
not  exist.  Within  the  literature  definitions  of  discontinuous 
behavior  such  as  “linear  kinematic  features”  (LKFs)  [27], 
[35],  “piece-wise  rigid  motion”  [35]  and  “aggregate  scale” 
[23],  [37]  are  beginning  to  emerge.  The  fundamental  prob¬ 
lem  in  defining  these  features  is  that  their  composition  tran¬ 
scends  any  one  scale,  with  several  scales  embedded  upon 
one  another. 

Fundamentally,  aggregate-scale  sea  ice  is  a  piece-wise 
continuous  material  made  up  of  a  collection  of  aggregated 
plates  subject  to  the  basic  principles  of  continuum  mechan¬ 
ics  with  limitations  of  those  basic  principles  in  close  prox¬ 
imity  to  discontinuities.  This  definition  places  the  motion 
of  sea  ice  under  the  larger  category  of  non-rigid  motion 
analysis  [26].  The  interdisciplinary  field  of  non-rigid  mo¬ 
tion  involves  the  understanding  of  three  basic  types  of  ma¬ 
terials,  namely,  continuous,  piece-wise  continuous,  and  dis¬ 
crete  particle  motion.  Depending  on  the  scale,  each  of  these 
descriptions  applies  to  sea  ice.  At  the  large  (basin)  scale, 
sea  ice  is  traditionally  regarded  as  a  non-rigid  continuum 
and  at  the  small  scale  it  is  regarded  as  a  collection  of  dis¬ 
crete  particles/floes. 

While  progress  has  been  made  in  understanding  these 
two  outer  scales  of  sea  ice,  the  intermediate  scale  of  piece- 
wise  continuous  sea-ice  motion  is  particularly  difficult  be¬ 
cause  it  goes  beyond  the  limits  of  its  scale  with  very  narrow 
(i.e.,  small  scale)  crack-like  features  stretching  thousands 
of  kilometers  (i.e.,  basin  scale).  This  scale  is  one  which  re¬ 
searchers  have  only  recently  been  able  to  explore  thanks  to 
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the  availability  of  high-spatial  resolution,  all-weather  Syn¬ 
thetic  Aperture  Radar  (SAR)  images  [15],  [16],  [17],  [19], 
This  paper  delves  into  the  estimation  and  validation  of  sea- 
ice  motion  using  SAR  images  from  the  European  Space 
Agency’s  (ESA)  first  European  Remote  Sensing  (ERS-1) 
satellite,  specifically,  imagery  from  C-band  (5.7  cm)  mi¬ 
crowave  SAR.  Among  the  algorithms  for  motion  tracking 
in  the  Antarctic,  the  method  by  Drinkwater  and  Kottmeier 
[18]  is  one  of  the  more  successful  ones.  The  data  set  ob¬ 
tained  using  their  method  from  Ice  Station  Weddell  (ISW) 
1992  [16],  [15]  is  the  one  we  chose  to  validate  our  algo¬ 
rithm. 


2.  Relevant  Background 

Since  the  original  description  of  the  problem  by  Horn 
and  Schunck  in  the  seminal  work  [24],  optic  flow  estimation 
has  been  a  much  researched  topic.  Under  the  assumption  of 
extremely  small  temporal  resolution  the  optic  flow  equation 
is  considered  valid  and  many  techniques  have  been  devel¬ 
oped  to  estimate  the  flow  field  [25],  [3],  [5],  Robust  tech¬ 
niques  [8],  [6],  [36]  have  emerged  to  handle  the  large  noise 
and/or  the  failure  of  the  underlying  image  motion  model. 

Sequential  pairs  of  SAR  images  capture  two  important 
characteristics  of  sea-ice  motion,  global  translation  and  dif¬ 
ferential  non-rigid  dynamics.  Due  to  differences  in  satellite 
orbits  and  the  earth’s  rotation,  the  visual  capture  of  local 
sea-ice  dynamics  is  overshadowed  by  the  complexity  to  as¬ 
certain  high  magnitude  global  translation  (on  the  order  of 
about  200  pixels  on  100m  resolution  images  separated  in 
time  by  three  (3)  days).  Traditionally  this  problem  has  been 
addressed  using  a  hierarchical  framework  [11],  [4],  Specif¬ 
ically  within  the  field  of  sea  ice,  this  problem  is  currently 
addressed  using  a  variety  of  methods  including  cross  cor¬ 
relation  and  2D  wavelets  [21],  [28],  [29],  [20],  [16],  [15], 
[31],  [30], 

An  added  complexity  in  this  traditional  methodology 
comes  from  the  fact  that  high-spatial  resolution  SAR  data 
is  limited  by  low-temporal  resolution  due  to  polar  orbital 
constraints  (typically  3  days).  Under  the  influence  of  fast 
moving  storms,  significant  non-linear  changes  in  disconti¬ 
nuities  can  occur  at  temporal  scales  much  lesser  than  3  days 
and  sea  ice  can  deform  rapidly  resulting  in  large  changes 
in  the  orientation,  distribution,  and  size  of  continuous  and 
discontinuous  regions.  Thus  estimation  of  the  motion  field 
in  SAR  imagery  requires  an  added  perspective  when  com¬ 
pared  to  the  traditional  optical  flow  algorithms.  The  hy¬ 
pothesis  considered  here  is  that  motion  can  be  extracted  in 
a  hierarchical  fashion  with  coarse  resolution  levels  resolv¬ 
ing  large  global  translation  through  linear  models  and  finer 
resolution  levels  resolving  smaller  local  non-rigid  dynamics 
using  higher  order  parametric  motion  models  such  as  affine, 
quadratic,  etc. 


3.  Motion  Estimation  Technique 


In  this  paper,  we  use  phase  correlation  as  the  technique 
to  estimate  sea-ice  motion.  Extensive  studies  on  the  phase 
correlation  process  for  motion  estimation  have  been  pro¬ 
vided  in  [40],  [41],  Various  image  registration  techniques 
using  the  Fourier  Transform  exist  in  current  literature  [12], 
[38],  [32],  Alliney,  et  al.  [2]  estimate  object  motion  under 
the  constraint  that  the  background  does  not  undergo  signif¬ 
icant  changes.  Stone,  et.  al.,  [39],  and  Foroosh,  et.  al., 
[22]  provide  algorithms  to  perform  subpixel  image  regis¬ 
tration  using  the  phase  correlation  technique.  Phase  cor¬ 
relation  as  described  in  research  literature  is  derived  from 
the  Fourier  Shift  Theorem  which  is  a  specialization  of  the 
Affine  Fourier  Theorem  as  proposed  by  Bracewell  [10]. 

The  affine  theorem  describes  the  separation  of  the  dis¬ 
placement  vector  and  the  affine  coordinate  matrix  in  the 
frequency  domain.  When  handling  the  large  motion,  as  ob¬ 
tained  in  the  satellite  imagery,  the  global  translation  compo¬ 
nents  overshadow  the  differential  affine  motion  parameters 
and  thus  the  affine  theorem  has  an  immense  potential  in  es¬ 
timating  the  motion  since  the  translation  component  can  be 
estimated  independent  of  the  affine  parameters. 

The  translational  component  can  be  extracted  from 
within  the  correlation  equation  using  ’’whitening”  zero- 
phase  FIR  filters  H\  =  ^(u)!”1  and  H2  =  |G*(u)|-1 
[33]  leading  to  a  Dirac  delta  function  centered  at  the  trans¬ 
lation  parameters  as  in  Eq.  1 


9.1(au)G^u))  =  iS(x_d) 


//|  Hn 


(1) 


The  main  advantage  of  phase-based  techniques  are  its 
characteristic  insensitivity  to  correlated  and  frequency- 
dependent  noise.  With  the  availability  of  the  2-D  FFT  [13], 
the  calculations  can  be  performed  with  much  lower  com¬ 
putational  complexity.  The  disadvantage  with  these  tech¬ 
niques  is  that  they  are  applicable  only  under  well-defined 
transformations  and  thus  require  bolstering  from  other  tech¬ 
niques,  especially  in  the  cases  where  the  transformations  are 
arbitrary.  A  point  of  consideration  is  that  phase  correlation 
peak  may  reduce  in  height  under  image  rotation  and  scal¬ 
ing  [32],  For  sea-ice  SAR  imagery,  these  can  be  assumed 
negligible  since  the  global  translation  is  large  compared  to 
affine  parameters  and  thus  not  significantly  affected  by  the 
magnitude  reduction. 


3.1.  Processing  Method 


Global  translations  of  sea  ice  in  SAR  imagery  are  of 
order  100  to  200  pixels  depending  on  the  temporal  dis¬ 
tance  between  images.  The  traditional  method  of  “Normal¬ 
ized  Cross  Correlation”  requires  a  large  support  window  to 


capture  large  translation  but  these  large  windows  encom¬ 
pass  a  combination  of  various  motions  such  that  accuracy 
of  the  estimation  procedure  is  limited  to  a  spatial  resolu¬ 
tion  much  larger  than  the  image  resolution.  Additionally, 
SAR  is  an  active  microwave  system  with  an  emitting  sig¬ 
nal  in  the  range  of  centimeters  (e.g.,  C-band  has  5.7  cm 
wavelength).  Materials  like  sea  ice  with  crystalline  features 
the  same  size  as  the  emitted  signal  introduce  considerable 
backscatter  which,  in  turn,  alters  and  blends  the  return  sig¬ 
nal  between  pixels.  Phase  correlation  has  the  basic  property 
of  being  illumination  invariant  and  can  therefore  be  used  ro¬ 
bustly  to  estimate  the  large  motion  at  a  low  computational 
burden. 

As  a  preprocessing  step  prior  to  the  estimation,  the  ERS- 
1  images  were  modified  using  histogram  equalization  in 
tandem  with  mid-tone  modification  to  obtain  ’’visually  sig¬ 
nificant  regions”  [7],  Experiments  with  various  equaliza¬ 
tion  methods  indicated  that  the  estimated  motion  field  had 
the  smallest  error  variance  when  mid  tone  expansion  was 
applied. 

To  tackle  the  large  translational  motion,  the  motion  field 
was  estimated  using  a  multi-resolution  hierarchy  of  images 
obtained  by  decimating  the  image  signal  using  a  median  fil¬ 
ter.  The  median  filter  was  used  instead  of  a  Gaussian  filter 
because  of  the  poor  response  of  Gaussian  Kernels  in  regions 
with  heterogeneity,  which  constitute  the  differential  features 
such  as  leads,  families  of  slip  lines,  cracks,  and  ridges. 

Each  image  resolution  hierarchy  obtained  from  an  image 
pairs  is  used  to  estimate  global  motion  by  performing  phase 
correlation  at  each  level  of  the  hierarchy  and  then  passing 
that  coarse  estimate  on  as  an  initial  estimate  to  the  next  finer 
level  of  the  pyramid.  The  global  motion  is  obtained  at  a  par¬ 
ticular  level  of  the  pyramid  by  binning  the  motion  vectors  at 
that  level  into  a  histogram.  The  histogram  bucket  containing 
the  maximum  number  of  candidates  is  taken  as  the  estimate 
at  that  level  of  the  hierarchy  and  submitted  to  the  next  level 
as  the  initial  guess.  This  process  is  repeated  down  to  the 
finest  resolution  image  in  the  pyramid  hierarchy  where  the 
final  estimated  solution  is  chosen  as  the  global  translation. 
To  improve  robustness  of  the  method,  potential  median  can¬ 
didates  at  each  level  in  the  pyramid  hierarchy  are  selected 
and  passed  through  a  “Lorentzian  estimator”  [9]  to  deter¬ 
mine  the  best  possible  estimate. 

Due  to  the  periodic  nature  of  the  Discrete  Fourier  Trans¬ 
form,  the  maximum  measurable  estimate  using  the  Fourier 
Transform  of  a  signal  within  a  window  of  size  W  is  W/2. 
Thus,  to  capture  translations  of  magnitude  (u,v),  the  win¬ 
dow  size  should  be  at  least  2  *  max  (it,  v).  For  the  ERS-1 
experimentation,  the  block  size  was  taken  as  32x32  and  the 
window  size  was  taken  as  128  x  128.  The  sizes  of  the  win¬ 
dow  and  the  block  are  maintained  constant  throughout  the 
entire  pyramid  hierarchy.  This  provided  a  method  of  ob¬ 
taining  accurate  estimates  at  the  coarsest  level  of  the  pyra¬ 


mid  which  amplified  to  large  translations  at  the  finer  scales 
of  the  image  pyramid.  Experience  with  ERS-1  SAR  im¬ 
ages  shows  that  a  three  level  pyramid  ranging  in  resolutions 
from  1536  x  1536  to  384  x  384  is  sufficient  to  capture  the 
observed  large  translation. 

Having  obtained  the  global  motion  compensated  image, 
a  local  piece-wise  linear  motion  model  is  applied  using 
phase  correlation.  This  local  motion  model  then  provides 
two  important  pieces  of  information  to  the  affine  paramet¬ 
ric  model:  (1)  an  initial  estimate  for  the  affine  motion  and 
(2)  a  correlation  map. 

For  this  experiment,  the  affine  model  chosen  is  a  least- 
square  regression  technique.  It  is  applied  making  use  of  a 
known  limitation  of  correlation  methods,  namely,  their  in¬ 
ability  to  resolve  regions  of  discontinuity  [41].  Low  values 
of  correlation  coefficient  on  the  correlation  map  are  thus 
used  to  isolate  discontinuities  so  that  the  motion  model  can 
track  the  deformations  in  the  vicinity  of  these  long  narrow 
crack-like  features  but  avoid  computations  right  at  the  dis¬ 
continuities. 

3.2.  Validation  Data  Set 

The  European  Space  Agency’s  ERS  -  1  and  ERS  -  2  C- 
band  (5.3  GHz,  5.7  cm)  Active  Microwave  Instrument  are 
currently  being  used  to  generate  weather  independent  (day 
or  night),  frequent  repeat,  high  resolution  (10  m)  100  km 
swath  radar  images  globally  and,  in  this  particular  case,  for 
the  Weddell  Sea  in  the  Southern  Ocean  around  Antarctica. 
The  5  month  Ice  Station  Weddell  1992  (ISW)  is  the  only 
winter  field  experiment  to  date  conducted  in  the  Western 
Weddell  Sea.  The  orbit  phasing  of  the  ERS  -  1  was  fixed  in 
the  3-day  exact  repeating  orbit  called  the  “ice-phase”  orbit 
providing  uninterrupted  SAR  imagery  of  spatial  coverage  of 
100  x  100  km2  during  the  entire  duration  of  the  experiment. 

The  SAR  images  obtained  from  ERS  -1  are  projected 
onto  the  polar  stereographic  Special  Sensor  Microwave  Im¬ 
ager  (SSM/I)  grid.  The  images  are  then  block  averaged  to 
100  m  resolution  (8x8  pixel  block  averaging)  to  speckle  fil¬ 
ter  the  images  to  minimize  intensity  errors  to  ±  1  dB  (>90% 
confidence  interval)  leading  to  images  with  dimensions  of 
1536  pixels  in  the  horizontal  and  vertical  direction. 

Motion  vectors  for  each  100x100  km2  SAR  image  are 
resolved  using  a  nested  cross-correlation  procedure  [16], 
[15]  to  characterize  5x5  km2  spatial  patterns.  A  total  of  1 1 
such  image  pairs  exist  from  this  processing  with  an  RMSE 
less  than  0.5  cm/s  for  3-day  tracking  compared  against  six 
buoy  measurements  from  the  ISW  [16]. 

4.  Results  and  Analysis 

The  global  estimates  obtained  using  our  processing 
method  are  compared  against  the  ISW  results  described  in 


Section  3.2.  The  global  motion  field  used  in  the  comparison 
is  obtained  by  converting  the  motion  field  from  pixel  coor¬ 
dinates  to  the  Special  Sensor  Microwave  Imager  (SSM/I) 
grid  coordinates  and  plotting  them  in  tandem  with  the  ISW 
motion  vectors. 

To  maximize  the  throughput  of  the  FFT  modules  the 
phase  correlation  is  performed  on  a  window  size  with  pow¬ 
ers  of  2.  Depending  on  the  spatial  resolution  required  for 
the  interpolation  of  the  motion  field,  the  block  size  can  be 
adjusted  at  8x8,  16x  16  or  32x32  leading  to  an  output  mo¬ 
tion  field  at  0.8  km,  1.6  km  or  3.2  km  resolution. 


1300 


1280 


1260 


1240 


•g.  1220 

£ 

9) 

g  1200 
w 

5  1180 


1160 


1140 


1120 


-1710  -1690  -1670  -1650  -1630  -1610  -1590  -1570  -1550  -1530  -1510 
Distancefkm] 


Magnitude  (kms) 


JPL  Vectors 


Direction  (degs) 


Figure  1 .  Comparison  of  motion  vectors  from 
the  ERS-1  image  pair  at  orbit  3402  frame  5103 
and  orbit  3412  frame  5693  shown  with  our 
estimated  vectors  (red)  and  the  ISW  vectors 
(blue).  Scatter  plot  comparisons  of  the  mag¬ 
nitude  and  phase  between  these  two  data  sets 
shown  in  accompanying  panels. 


The  subsequent  results  in  Figure  1  show  the  accuracy  of 
the  estimates  obtained  by  the  phase  correlation  technique 


in  comparison  with  the  ISW  vectors.  The  estimated  mo¬ 
tion  field  in  the  images  below  have  been  computed  using 
a  32x32  window  and  then  overlaid,  using  linear  interpola¬ 
tion,  on  a  5  km  SSM/I  grid  for  comparison  with  the  ISW 
vectors.  Scatter  plots  show  the  magnitude  and  phase  vari¬ 
ation  between  our  estimates  and  the  ISW  vectors  for  the 
spatially  collocated  positions. 

Figure  1  is  an  illustrative  example  of  the  comparison  be¬ 
tween  the  ISW  ground  truth  vectors  (in  blue)  and  the  esti¬ 
mated  vectors  (in  red).  The  bottom  panel  in  Figure  1  shows 
scatter  plots  of  the  magnitude  and  directional  components 
of  the  estimated  vectors  with  the  ISW  vectors. 

Two  statistical  measures  of  the  similarity,  the  Root  Mean 
Square  Error  (Eq.  2)  and  the  index  of  agreement  (Eq.  3) 
have  been  computed. 


RMSE  = 


N 


Jc  £(p*  -  °fc)2 

k= 1 


(2) 


.  =  ,  _  Ek^iMPk  -pfc)7 

7  i:LiM\Pk-o\  +  \ok-o\p 

where  pt  are  the  predicted  samples,  Ofc  are  the  observed 
ground  truth  vectors,  are  the  weight  functions  which 
are  assumed  uniform  for  this  study  and  o  is  the  mean  of 
the  observed  data.  7  is  the  order  of  index  and  according  to 
Willmott  [43],  7  =  1  is  most  robust  for  comparing  results 
because  of  its  linear  approach  to  a  perfect  match. 

Analysis  of  the  image  pairs  using  higher  resolution 
analysis  windows  (Figure  2)  reveals  precise  demarcations 
corresponding  to  the  regions  in  the  ice-flow  undergoing  sig¬ 
nificant  amounts  of  non-rigid  dynamics.  Based  on  the  vari¬ 
ance  of  the  magnitude  and  direction  of  the  motion  field,  the 
cluster  map  as  shown  in  Figure  2  was  created  using  a  quad¬ 
tree  model. 

This  result  reveals  the  usefulness  of  local  higher  order 
motion  in  the  vicinity  of  the  regions  of  discontinuity.  Given 
the  global  motion  compensated  images,  local  non-rigid  dy¬ 
namics  was  first  extracted  using  the  simplest  model  of  the 
local  motion.  Under  the  assumption  that  magnitude  of  the 
differential  motion  is  small,  a  piecewise  linear  approxima¬ 
tion  of  the  non  rigid  motion  using  phase  correlation  was 
applied.  As  seen  in  Figure  3,  the  estimated  local  motion 
vectors  are  overlaid  upon  the  correlation  map  to  verify  the 
goodness-of-fit  as  are  the  ISW  vectors  in  validation.  Exper¬ 
imentally,  a  threshold  correlation  of  <  0.2  is  found  to  be  the 
most  appropriate  for  demarcating  between  continuous  and 
discontinuous  regimes.  This  low  correlation  corresponds  to 
regions  where  the  magnitude  of  deformation  far  exceeds  the 
piece-wise  approximation. 

The  dynamic  region  analyzed  and  described  earlier  in 
(Figure  2)  is  shown  again  in  Figure  4  in  the  same  format 
as  Figure  3  for  comparison.  This  dynamic  region  contains 
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Figure  2.  Comparison  between  the  estimated  vectors  and  ISW  vectors  from  the  SAR  image  pair  at 
orbit  3197  frame  5693  and  orbit  3230  frame  5103  at  a  dynamic  zone. 


lots  of  pixels  without  vector  solutions  because  these  pixels 
did  not  meet  our  correlation  threshold.  Such  areas  require 
a  higher  order  model  than  our  method  currently  includes. 
We  note  however,  that  the  local  motion  is  accurate  in  the 
high  correlation  regions  as  validated  through  comparison 
with  the  ISW  vectors. 

Due  to  the  projection  of  the  non  linear  components  of  the 
higher  order  motion,  the  local  motion  so  estimated  contains 
false  discontinuities  (i.e.  false  positives).  These  discontinu¬ 
ities  occur  due  to  abrupt  changes  in  the  frequency  compo¬ 
nents  between  the  two  regions  under  consideration  causing 
variations  in  the  estimated  vector  field.  To  reduce  the  effect 
of  these  false  discontinuities,  sub-pixel  motion  interpola¬ 
tion  was  carried  out  using  a  cubic  spline  within  a  window 
around  the  result  of  the  local  phase  correlation.  This  proce¬ 


dure  reduced  the  bands  of  discontinuities  within  the  motion 
field  and  also  estimated  the  motion  field  at  a  finer  resolu¬ 
tion  around  the  peak  of  the  local  estimate  but  at  the  cost  of 
increased  computational  complexity. 

Figure  5  is  the  combined  scatter  map  over  all  the  im¬ 
age  pairs  used  in  this  analysis.  As  can  be  seen,  the  ma¬ 
jority  of  data  points  (estimated  vectors)  cluster  along  the 
one-one  line  with  very  small  dispersion.  ncornpiete  is  the  to¬ 
tal  number  of  motion  vectors  and  ngood  are  the  points  that 
are  within  400  m  (4  pixel)  of  the  ISW  vectors.  This  devi¬ 
ation  is  at  the  threshold  of  SAR  geolocation  accuracy  [15] 
and  therefore  a  scatter  based  on  instrument  accuracy  rather 
than  the  technique  applied,  rideformation  are  all  the  points 
in  the  correlation  maps  having  low  correlation.  Thus  72.8% 
of  the  estimated  motion  vectors  are  accurate  with  respect 
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Figure  3.  Relative  motion  field  between  SAR  image  pair  of  orbit  3412  frame  5693  and  orbit  3455  frame 
5693  is  shown.  Grey-scale  map  plotted  with  lightest  to  darkest  denoting  high  to  low  correlation 
coefficient  with  darkest  defining  discontinuity  regions  (without  motion  vectors)  in  the  high  resolution 
product  (0.8  km  resolution  thin  black  vectors).  Interpolated  motion  vectors  at  5  km  resolution  (thick 
black)  further  superimposed  with  coincident  ISW  vectors  (thick  white)  which  hide  nearly  all  of  each 
thick  black  vector  in  validation. 
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Figure  4.  Relative  motion  for  SAR  image  pair  with  orbit  3197  frame  5693  and  orbit  3230  frame  5103 
showing  the  highly  dynamic  region  presented  in  Figure  2  earlier. 
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Figure  5.  Scatter  plot  of  the  estimated  vec¬ 
tors  v/s  ISW  vectors  for  the  entire  sequence  of 
images  during  ISW.  Green  indicate  the  good 
points,  blue  are  the  deforming  points  and  red 
are  points  we  can  not  explain. 


to  ISW  vectors.  Including  flagged  points  undergoing  de¬ 
formation  we  can  account  for  93.5%  of  all  the  data  points. 
The  remaining  6.5%  of  data  are  displacement  results  greater 
than  400m  which  we  can  not  account  for  using  a  correlation 
map.  These  points  mainly  correspond  to  regions  where  the 
gradient  of  the  velocity  is  high  and,  in  principle,  requires 
a  higher-order  motion  model  to  localize  the  position  of  the 
discontinuity  accurately.  Also  indicated  in  the  scatter  plot 
are  the  RMSE  (Eq.  2)  and  index  of  agreement  (Eq.  3)  for 
the  complete  data  set  and  the  “good”  data  set  indicating  a 
high  accuracy  between  the  ISW  vectors  and  our  estimated 
motion  vectors. 

To  improve  upon  the  estimated  parameters,  an  affine  mo¬ 
tion  is  imposed  onto  the  local  piece-wise  linear  motion. 
This  is  performed  specifically  in  the  regions  of  discontinuity 


as  obtained  from  the  correlation  map.  The  affine  parameter 
estimation  chosen  is  a  least-squares  fit  to  the  locally  trans¬ 
lated  patch.  The  Peak  Signal  to  Noise  Ratio  (PSNR)  values, 
calculated  using  Eq.  4,  are  provided  in  Table  1. 

255 

PSNR  =  20  *logw[ - ]  (4) 
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where  <r  is  the  standard  deviation.  Image i  and  Image 2  is 
the  image  pair  under  consideration.  From  the  table  it  is  ev¬ 
ident  that  the  PSNR  value  between  Image 2  and  the  mo¬ 
tion  compensated  images  increases  with  increasing  motion 
model  complexity. 


5.  Conclusion 


Localizing  and  parameterizing  discontinuities  is  an  im¬ 
portant  landmark  for  sea-ice  research.  Due  to  the  low  tem¬ 
poral  resolution  of  the  SAR  images,  this  type  of  information 
is  paramount  input  into  numerical  models  for  interpretation 
of  events  between  SAR  scenes.  This  paper  describes  a  hi¬ 
erarchical  method  to  estimate  the  parameters  of  a  discon¬ 
tinuous  motion  field  such  as  the  ones  captured  by  the  ERS- 
1  SAR  imagery  in  three  stages.  The  hypothesis  of  using 
a  finer  sieve  to  filter  out  coarse  motion  models  iteratively 
seems  to  be  validated  with  the  observed  improvement  in  the 
PSNR  values  of  the  motion  compensated  images. 

The  comparative  results  between  the  estimated  vectors 
and  the  ISW  ground  truth  vectors  indicate  that  hierarchical 
phase  correlation  is  ideal  for  the  estimation  of  the  global 
motion.  This  is  mainly  due  to  the  inherent  robustness  of 
phase  correlation  to  illumination  variation  and  its  ability 
to  capture  large  translational  components  without  getting 
caught  in  local  minima.  Additionally,  the  usage  of  the  Fast 
Fourier  Transform  in  the  calculation  of  the  phase  correlation 
term  makes  the  computation  significantly  faster  then  optic 
flow  methods. 

The  second  stage  involves  isolation  of  discontinuous  re¬ 
gions  using  local  piecewise  linear  model  followed  by  the 
third  stage  of  estimating  affine  parameters  along  regions  of 
discontinuities.  Under  the  assumption  that  the  net  motion 
is  actually  composed  of  a  large  global  motion  component 
and  small  local  deformations,  this  method  captures  the  large 
global  motion  component  and  the  local  deformation  using 
an  affine  motion  model. 

As  a  subsequent  stage  to  the  current  research,  the  esti¬ 
mate  of  the  local  deformation  could  be  improved  using  ro¬ 
bust  parameter  estimation  techniques  and  also  inclusion  of  a 
motion  model  such  as  the  quadratic.  Another  possible  track 
for  future  research  would  be  a  feature-based  approach  in 
tandem  with  global  motion  estimation  in  order  to  improve 
the  overall  robustness  of  the  estimated  global  motion. 
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Table  1.  Comparison  of  PSNR  between 
motion  compensated  images  and  lmage2, 
PSNRg  -  Global,  PSNRLL  -  Local  Linear, 
PSNRLyl  -  Local  Affine 


Image  i 

Image2 

PSNRg 

psnrll 

psnrla 

2982-5693 

3025-5693 

13.109 

13.934 

14.116 

3058-5103 

3068-5693 

13.915 

14.834 

15.006 

3068-5693 

3111-5693 

14.865 

16.135 

16.477 

3111-5693 

3144-5103 

11.419 

11.611 

11.649 

3144-5103 

3154-5693 

14.128 

15.268 

15.414 

3154-5693 

3197-5693 

11.351 

11.395 

11.379 

3197-5693 

3230-5103 

15.715 

18.538 

19.080 

3230-5103 

3283-5693 

14.213 

15.645 

15.947 

3402-5103 

3412-5693 

12.489 

13.944 

14.078 

3412-5693 

3455-5693 

12.850 

13.645 

13.689 

3412-5713 

3455-5713 

13.446 

14.524 

14.699 
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